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Abstract 

When particles suspended in a fluid are driven through a regular lattice of cylindrical obstacles, 
the particle motion is usually not simply in the direction of the force, and in the high Peclet number 
limit particle trajectories tend to lock along certain lattice directions. By means of molecular 
dynamics simulations we show that this effect persists in the presence of molecular diffusion for 
nanoparticle flows, provided the Peclet number is not too small. We examine the effects of varying 
particle and obstacle size, the method of forcing, solid roughness, and particle concentration. While 
we observe trajectory locking in all cases, the degree of locking varies with particle size and these 
flows may have application as a separation technique. 
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I. INTRODUCTION 



Macroscopic particles suspended in a fluid and driven through a periodic lattice of solid 
obstacles tend to move in (commensurate) periodic trajectories. In general, these periodic 
trajectories are not aligned with the driving force. In fact, the average trajectory an- 
gle typically remains locked to a given lattice direction over finite intervals of the forcing 
angle (directional locking). The locking direction a (= arctan(g/p), where (p,q) is a lattice 
direction) and its relation to the forcing angle 9 depends on the properties of the particles. 
Therefore, particles of different species driven by the same external force could migrate at 
different angles with respect to the underlying lattice of solid obstacles. This behavior raises 
the possibility of using periodic obstacle lattices to provide deterministic and continuous 
separation of suspended particles. In fact, we have experimentally shown the possibility to 
separate sedimenting particles by size as they move through an ordered array of cylindrical 
obstacles created with LEGO® pegs on a LEGO® board. [2J 

Experiments on colloidal particles moving through periodic arrays have shown analo- 
gous directional-locking behavior in the near-deterministic limit of large Peclet numbers, jj] 
Moreover, a separation technique based on the deterministic motion of suspended particles 
in these periodic devices has been successfully used to fractionate colloidal particles of dif- 
ferent size, [3|, isolate blood plasma from other blood components separate and label 
cells [5j as well as to isolate selected cell types for tissue engineering applications. [6| Similar 
locking behavior is observed in the transport of suspended colloidal particles through spa- 
tially periodic, smooth, two-dimensional energy landscapes, arising from an array of optical 
tweezers, pj Introducing an optical lattice in a microfluidic channel has also been shown to 
selectively induce locked-in trajectories that deflect specific particle species in a direction 
different from the average flow, thus providing a means for lateral fractionation based on 
the dielectric properties of the particles. 0, 9] An analogous optical fractionation method is 
based on the presence of locked-trajectories depending on particle size. [10, 

Simulations which faithfully incorporate all of the relevant experimental features of col- 
loidal particles - finite particle and obstacle size, accurate hydrodynamics and Brownian 
motion - are not as yet available. The motion of tracer particles moving through an array 
of solid obstacles was recently investigated using a Fokker-Planck equation approach, and 
directional locking was observed at large Peclet numbers. 12[] However, these simulations 
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do not account for particle-solid hydrodynamic interactions or other particle size effects. 
The case of macroscopic particles was recently investigated by means of Stokesian dynamics 
simulations, showing that directional locking occurs in the presence of irreversible particle- 
obstacle interactions.pj These irreversible forces, on the other hand, are added ad hoc, 
typically to prevent overlap in the Stokesian dynamics simulations. In addition, these simu- 
lations correspond to the limit of zero inertia and deterministic motion. Directional locking 
(or phase-locking) has also been observed in simulations of a number of related transport 

n 

phenomena, inclu ding driven vortex lattices, [13J the motion of overdamped particles in pe- 
riodic force fields, 141 particles driven through a colloidal lattice, 151 particles moving on 



crystalline surfaces, 



161 ] and periodic energy landscapes in general. 17] Most of these sim- 



ulations, however, investigate the motion in an external field and not the interaction of 
suspended particles with solid obstacles. 

In this paper we use molecular dynamics simulations to show that directional locking 
occurs in the transport of nanoparticles, provided the Peclet number is not too small. 
These simulations allow us to examine the complex interplay between Brownian motion 
and different driving fields, together with hydrodynamic interactions (with possible inertia 
contributions), and other finite size effects such as the atomic-scale roughness of the parti- 
cles. In addition, molecular dynamics simulations also allow us to investigate the effect that 
particle-particle interactions have on the separation behavior at different particle concentra- 
tions. Moreover, at low concentrations and for certain orientations of the driving force, we 
show that particles of different size move in different directions. 



II. METHODS 



The calculations are based on standard molecular dynamics techniques [18 
interacting via Lennard- Jones (LJ) potentials: 

-12 



4e 



for atoms 



(1) 



cut off at r c = 2.5a. The coefficient cy will be used to vary the interaction according to 
the atomic species i and j involved. We have a three component system consisting of a 
monatomic suspending liquid (L), mobile nearly-spherical particles (P) composed of atoms 
rigidly held together, and cylindrical obstacles (O) composed of atoms fixed in place. The 
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FIG. 1: (Color online) Snapshot of the atomic positions for a simulation involving a 4 x 4 array 
of cylindrical obstacles of radius 2 (green) and a single mobile sphere (blue), also of radius 2, 
suspended in a monatomic solvent (red dots). Left: x-y plane in which the force is applied; right: 
x-z plane. 



interaction coefficients are chosen to have a "standard" liquid, en = 1, which completely 
wets both solids, clo = clp = 1, and only a short-distance repulsion between the particles 
and the obstacles, cop = 0. We express all dimensional quantities in terms of the energy 
scale e, the core size a, and the natural time scale r = o^m/e) 1 / 2 , where m is the common 
mass of all atoms. Typical values are e ~ 10~ 21 J, a ~ 0.3 nm and r ~ 2 ps. 

The initial atomic positions form a cubic lattice of density p = 0.8<j -3 . The particle is 
denned by selecting all atoms within a prescribed radius of some initial center, while the 
obstacles are centered on a square lattice in the x-y plane, and consist of all atoms whose 
(x, y) coordinates lie within a second prescribed radius of one of the centers. An example 
of the atomic positions after equilibration is given in Fig. [T], showing two orthogonal views 
of a 4 x 4 lattice of cylindrical obstacles of radius 2, along with a particle of radius 2. The 
simulated region is a cube 43.09<r on a side, containing 56,288 solvent atoms, 16 obstacles 
with 480 atoms each, and the mobile sphere contains 32 atoms. Because they are composed 
of sets of distinct atoms the particle and obstacle automatically have rough surfaces, and we 
discuss the effects of varying the degree of roughness below. Periodic boundary conditions 
are applied in all three directions; the spacing between the obstacles must be commensurate 



4 



with the size of the simulation box so that the obstacle lattice has the same periodicity, and 
in the case shown in the figure the ratio is 1/4. Although the particle could fully sample the 
periodic geometry in a smaller simulation involving just one unit cell of the obstacle lattice, 
the larger system provides a larger separation between the particle and its periodic images 
and the particle's motion is correspondingly less influenced by the "wake" of the images. 

The motion of the atoms follows the Newtonian equations appropriate to each species. 
Each fluid atom accelerates according to the force exerted by all other atoms within inter- 
action range. The particle moves rigidly, so that its center of mass translates under the sum 
of the LJ forces applied to its individual atoms and rotates in response to the sum of the 
torques on the individual atoms about the center of mass. The orientation of the particle 
and its rotational motion is conveniently described by the Euler equations expressed in terms 
of quaternion variables. [19( The atoms in the obstacles are simply fixed in place. The par- 
ticle is driven through the obstacle array by adding an external force to its LJ interactions 
in a prescribed direction in the x-y plane, or alternatively such a force is applied to each 
suspending fluid atom to produce a pressure driven flow in the solvent in which the particle 
is advected. 

During the simulation the temperature of the solvent liquid is maintained at T = l.Oe/fcg 
by a Nose-Hoover thermostat. In principle, the atoms in the particle and obstacle should also 
have thermal fluctuations about their rigid-lattice positions, but the displacements in a solid 
are small and this complication has been omitted in most cases. Allowing thermal motion 
of the obstacle atoms presents little computational difficulty, and can be accomplished by 
using a stiff linear spring to tether these atoms to their initial lattice sites and allowing them 
to move in response to the LJ interactions with other nearby atoms. We have added such 
thermal motion to the obstacles in a few test cases and, outside of an interesting special 
situation discussed below, we do not observe any significant change in the particle motion. 
Adding thermal motion to the atoms of the mobile particle is more difficult because the 
particle's moment of inertia would vary in time and complicate the solution of the Euler 
equations. 

In these calculations the directed motion of the particle is modified by the Brownian 
motion resulting from its interaction with the suspending fluid atoms, and the competition 
may be quantified by the Peclet number. Normally one defines Pe = Ua/D m where U, a and 
D m are the average velocity, radius and molecular diffusivity of the particle, respectively, but 
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this form is awkward here because the velocity is a result of the simulation and measuring the 
diffusivity in this geometry would require extensive additional simulations. (The motion in z 
is diffusive, but the diffusion coefficient is modified from that in pure fluid by the interaction 
with the rough obstacles.) However, a simpler characterization is possible if we assume that 
a simple linear drag law applies, U = 7F, and that the drag coefficient satisfies the Stokes- 
Einstein relation D m = ksT/^y, in which case Pe = Fa/keT. In fact, as we shall see below, 
7 depends on the orientation of F with respect to the obstacle lattice so this argument is 
only approximate, and the latter form for Pe is just a convenient parametrization. 

III. RESULTS 

We begin with simulations of the system shown in Fig. [TJ applying a force 

F = (F cos 9, F sin 9, 0) (2) 

to the particle. The position of the particle's center of mass is recorded at lr intervals, and 
the subsequent analysis is based on the x and y components of these trajectories. Because of 
the symmetry of the square obstacle lattice, it suffices to consider forcing angles 9 between 
and 45°. 

A. Local behavior 

The effects of molecular diffusion are illustrated in Fig. [2] which shows a fine-scale view of 
the initial particle motion at low and high Fq. As expected, the motion between obstacles 
tends toward a straight line, and the motion around obstacles follows their shape more 
closely at larger Fo] note that numerically Pe = 2Fo here. At larger scales, the shape of the 
trajectories is quite sensitive to locking phenomena. In Fig. [3] we display sample trajectories 
at Fo = lOOe/a for forcing angles 9 = 36, 30, 26 and 22°, separately for intermediate and long 
time intervals of about 500 and 5000r, respectively. When the force direction is close to a 
locking angle of 45° (not shown), the particle follows a straight path with detours around the 
obstacles, and when the forcing angle 9 is near zero (not shown) the particle bounces back 
and forth between two rows of obstacles. These two limiting behaviors persist for angles 
somewhat larger than 0° and somewhat smaller than 45°, respectively. At intermediate 
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FIG. 2: (Color online) Early stages of the trajectories of a particle forced at 9 = 48°, shown as a 
sequence of circles (red) moving through cylindrical obstacles (green). Left: Fq = 5 until 2250r; 
Right: F = 200 until 25r. 

angles in the middle of this interval the motion is rather irregular: segments of variable 
length along the locking directions at 45, 26.6 and 0°. These directions correspond to 
displacements which would geometrically intersect the obstacle lattice with 8y : 5x of 1:1, 
1:2 and Loo; the other lattice possibilities such as 1:3 or 2:3 do not seem to appear in the 
trajectories to any significant degree, even when the force is applied at just those values. 

At this point we note an important negative result: locking in three-dimensional flows 
requires cylindrical obstacles which fully obstruct the motion in the third dimension, whereas 
a cubic array of localized obstacles has no systematic effect on the motion. The reason is that 
the deflection of the sphere around an obstacle after a close approach should be independent 
of z in order for the process to repeat coherently at subsequent obstacles. A cylindrical 
obstacle has this feature, but if the obstacle is spherical the impact parameter depends on 
z and will vary as the particle diffuses in that direction. The direct evidence is given in 
Fig. H] which shows the x-y trajectories of a mobile sphere of radius 2a moving through a 
cubic array of spheres of the same radius spaced by 10.77a. A force given by Eqj2]is applied 
with F = 50 e/a at forcing angles 9 = (22.5, 24.75, 27, 29.25, 31.5) and the inclination of the 
resulting trajectories is essentially the same: a = (22.27,24.51,27.08,29.19,31.52). In this 
case a single realization at each angle suffices. 
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FIG. 3: (Color online) Initial and global trajectories of a particle forced at Fq = 100e/a for angles 

8 = 36 (red), 30 (blue), 26 (cyan) and 22° (magenta), top to bottom. 

2000i I 

1800- 
1500- 
1400- 




FIG. 4: (Color online) Trajectories of a sphere forced through a three-dimensional cubic ar- 
ray of spherical obstacles, for forcing angles 6 = (22.5,24.75,27,29.25,31.5) bottom to top 
(red, blue, cyan, magenta, green). No locking occurs in this case. 

B. Locking behavior 

The locking behavior of the particle motion is determined by the relation between the 
forcing angle 9 and the trajectory angle a, defined as the angle between the end-to-end 
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vector and the x-axis. This angle could depend on the length of the trajectory, but in these 
simulations we observe little sensitivity once the particle covers distances beyond a few 
hundred a. The results presented here are based on simulations over a 5000r time interval 
for flows with Pe > 100, which cover at least this distance. For each case we average over 
five realizations, obtained by applying a force in the (neutral) z-direction for different short 
intervals, to produce a new initial position for the particle. The resulting standard deviation 
in the mean trajectory angle is one degree or less, and similarly the relative fluctuations in 
the other quantities discussed below is small. We begin with simulations in the system 
discussed above, using forces of magnitude F = 50, 100 or 200e/cr oriented at angles 14 
through 42°. As an alternative to applying a fixed force to the particle which pulls it through 
the solvent, in a separate series of simulations we apply a pressure gradient to the fluid which 
drives a flow and advects the particle. In the latter case, a force g is applied to each fluid 
atom at the same set of directions, with magnitude g = 0.1, 0.25 or 0.5cr/V 2 . 

The resulting variation of the trajectory angle is shown in Fig. [5], where obvious locking 
behavior is seen: all of the curves depart significantly from a straight line at 45°. For smaller 
values of 9 the trajectory angle is approximately zero in all cases. As expected, the locking 
behavior is more pronounced at higher forcing values, corresponding to higher Pe. A more 
surprising result is that a locking plateau - a range of forcing angles producing the same 
trajectory angle - appears only at a subset of the possibilities which occur in the ballistic 



and high Pe limits, jl, 12j at zero, 45 and possibly 26.6°. The latter near-plateau is most 
evident when the solvent is forced, at the two higher g values. Presumably the difference 
arises from molecular diffusion, even for Peclet numbers of several hundred which appear at 
the highest forcing values here. 

When higher values of force are applied to the particle a complication arises: the particle 
can become jammed into an obstacle and remain there for a variable length of time. The 
reason for jamming is that both particle and obstacle are rigid solid lattices with (atomic) 
corners, and if a corner of the mobile particle enters the gap between two layers of the obstacle 
and the driving force exceeds the typical LJ short-distance repulsive force, an equilibrium 
may be established. A simple (physically realistic) remedy is to introduce thermal motion 
for the atoms in the two solids, which alters the force balance and mobilizes the particle. 
This technique was used in a few test situations at F = 200, with no significant change 
in the results, and could be used to reach higher Pe values, but we have not pursued this 
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FIG. 5: (Color online) Trajectories angles a vs. forcing angle 9 when the particle is driven either 
by a constant force applied to the particle of magnitude Fq = 50 (green,*), 100 (cyan,o) or 200e/cr 
(magenta,x), or else by an acceleration g = 0.1 (red,-), 0.25 (blue,+) or 0.5cr/r 2 (orange,-^) applied 
to the solvent atoms. 

extension systematically. 

In addition to the angle we measure the average velocity V in each trajectory, defined as 
the end-to-end distance in the x-y plane divided by the elapsed time, and an "excess length" 
£, defined as the (approximate) path length divided by the end-to-end distance. The path 
length is approximated by dividing the trajectory into straight segment of duration It and 
summing the (two-dimensional) segment displacements, and the significance of £ is that it 
measures the deviation of a trajectory from a straight line and thereby the degree to which 
an individual trajectory is in fact locked at the angle a. The value of £ of course depends 
on the time interval of the segments used to measure it and is not unique, but rather is 
meant to provide a simple measure of "dithering" - deviations from a straight line along a 
trajectory. The result for the average velocity and excess length are presented in Fig. [6] The 
velocities are largest and the excess length smallest when the forcing direction is along the 
x-axis and the particle simply moves through one row of the lattice, and similar behavior 
is found when the forcing is close to the 45° locking direction. In the forced-particle cases, 
the velocity is smallest and the dithering is largest at intermediate forcing angles where the 
particle trajectory does not follow a single locking angle. The forced-solvent simulations 
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FIG. 6: (Color online) Average velocity (left) and excess length (right) for the trajectories in Fig. [5J 
plotted with the same colors and symbols. 

at g — 0.25 and 0.5<t/t 2 show some evidence for a plateau at the 26.6° locking angle, and 
correspondingly have local maxima in V and local minima in £. 

Next, we consider the effect of varying the geometry. We compare the behavior of particle 
of radius 2, 4 and 5.5<r, in a lattice of cylinders of the previous radius 2 but a larger spacing 
17.2a in order to accommodate the larger particles. A force Fq = 100 is applied to the 
particle in each case. In Fig. [7| we show the trajectory angle, average velocity and excess 
length as a function of forcing angle. In this situation, we again see that the trajectories do 
not simply follow the force direction, and at least for the two smaller particles the presence 
of three locking plateaus. In comparison to the previous case, the larger obstacle spacing 
leads to a 26.6° plateau for particles of radius 2, which suggest that a looser obstacle lattice, 
or more precisely a larger ratio of obstacle spacing to particle size, may improve locking. The 
overall velocity decreases with particle size, reflecting the increase in the hydrodynamic drag 
coefficient with radius, and the variation with angle shows similar trends as in the previous 
geometry, high velocities and low £ near the locking direction at 45° and near 26.6° for the 
two smaller particle sizes where a locking plateau is present. The radius 4 simulations, which 
extend into the 0° locking regime, show a velocity maximum and an excess length minimum 
there, consistent with the previous system, but the other particle cases have not covered 
this regime. 
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FIG. 7: (Color online) Effects of size. Trajectory angle (top), average velocity (left) and excess 
length (right) as a function of forcing angle for particles of radius 2 (red,+), 4 (green,*) and 5.5 
(cyan,o) in the larger-spacing obstacle lattice. 

We have also varied the particle and obstacle size so as to explore the effects of surface 
roughness. The scale of the roughness is controlled by the size of the atoms comprising the 
particle and obstacles, so if we double their radii we have reduced the relative roughness by 
a half. On the other hand, if at the same time we double the obstacle spacing we nominally 
have scaled the system up in size by a factor of two, and to the extent that these simulations 
are in the Stokes flow regime (the suspending fluid's viscosity is about 2m/ (err) so in these 
simulations the Reynolds number is 0(1)), the hydrodynamics is unchanged. In Fig. Owe 
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FIG. 8: (Color online) Effects of roughness. Trajectory angle (top), average velocity (left) and 
excess length (right) as a function of forcing angle for particles and obstacles of radius 2 spaced 
by 10.77 (green,*) at force Fo = 50, compared to a system doubled in all sizes at force Fq = 100 
(red,o). 

compared the results of doubling all sizes and applying a force of magnitude 100 to the 
original radius 2 case at a force of 50; doubling the particle size approximately doubles the 
hydrodynamic drag so as to provide a fair comparison. As seen in the figure, there is no 
significant change in the trajectory angle, and approximate agreement in the velocities and 
excess length. Theses results suggest that it is the absolute scale rather than the relative 
amount of roughness which controls the behavior. 
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FIG. 9: (Color online) Trajectory angles a as a function of the forcing direction 9. The solid 
lines correspond to Stokesian dynamics simulations with roughness parameters as indicated. The 
symbols correspond to the molecular simulations in Fig. [5]for Fq = 200e/c (solid) and g = 0.5u/r 2 
(open). 

C. Deterministic behavior in the Stokes limit 

In previous work we investigated the high Pe limit by means of Stokesian dynamics 
(SD) simulations, and showed the presence of directional locking. [lj] Stokesian dynamics 
simulations incorporate hydrodynamic interactions between the suspended particles and the 
solid obstacles in the limit of zero Reynolds number. In addition, a short-range repulsive 
force is typically introduced in SD simulations to qualitatively model the effect of non- 
hydrodynamic interactions, such as irreversible forces originating from surface roughness. 
In fact, the observed locking behavior was shown to depend on a small parameter e that 
determines the range of the repulsive force or the relative magnitude of the surface roughness. 

Here, we compare SD simulations with the behavior of nanoparticles at relatively large 
Peclet numbers. The obstacle lattice in the SD simulations is represented by an array of fixed 
spheres and the suspended particle moves in the plane of the array (note that in the case of 
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FIG. 10: (Color online) Excess length £ as a function of the forcing angle 9 calculated from Stokesian 
dynamics simulations for different values of the roughness parameter e as indicated. The transitions 
between locking directions are indicated in the figure, where (p, q) corresponds to motion at an 
average angle a = arctan(g/p). 

deterministic trajectories we do not have the problem discussed in section IIII Al regarding 
the absence of locking in the case of spherical obstacles given that the particle does not 
diffuse in the z direction). The particles forming the obstacle lattice as well as the moving 
sphere have the same radius a = 1 and the lattice spacing is 5a, analogous to the system 
used in the nanoparticles simulations. 

The locking behavior is presented in Fig. [9] for different values of the roughness param- 
eter, ranging from e = 5 x to e = 5 x 10~ 2 . Typically, particles with larger roughness 
exhibit fewer locking angles, as seen in the figure. We also compare the SD results with 
the trajectory angles obtained for nanoparticles shown in Fig. [5] for the largest values of 
the Peclet numbers. We observe that the migration angle of the nanoparticles is closer to 
the deterministic trajectories corresponding to the largest value of the roughness parame- 
ter. It is also clear that the gravity driven nanoparticles exhibit a plateau similar to that 
observed in the SD simulations for intermediate forcing angles. On the other hand, the 
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FIG. 11: Average velocity (top) and excess length (bottom) for the trajectories shown in Fig. [10] 
for the smallest value of the roughness parameter, e = 5 x 10~ 4 . 

force-driven nanoparticles seem to transition directly from a. = 0° to a = 45°. This behav- 
ior corresponds to deterministic trajectories with values of the roughness parameter larger 
than those considered here. 

The variation of the excess length measured in the SD simulations is also analogous to the 
one observed in the case of nanoparticles, but exhibits sharper transitions between locking 
angles. In Fig. [10] we present the excess length as a function of the driving direction for 
the different values of the roughness parameter considered in Fig. [9] Similar to what we 
observe in the transport of nanoparticles at large Pe, the excess length is smaller when the 
particles move either at a = 0° or a = 45°. At intermediate angles the excess length displays 
sharp jumps or local maxima at the forcing angles corresponding to transitions between two 
locking directions. Similar dependence is observed in the motion of nanoparticles presented 
in Figs. [6] and [7] In the cases in which there is a plateau at intermediate driving angles 
the excess length displays two clear maxima. On the other hand, in the cases in which 
a = 0° and 45° are the only locking directions, there is a single maximum in the excess 
length around the transition angle. 
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The deterministic nature of the SD simulations also allowed us to investigate the observed 
correlation between the average velocity and the excess length in the trajectory of the par- 
ticles. In Fig. [TT]we see that, in fact, they are highly correlated, with local maxima in the 
excess length corresponding to minima in the average velocity of the suspended particles. In 
fact, the slowest regions of a given trajectory occur when the particle moves around a solid 
obstacle with small particle-obstacle separations. On the other hand, the fastest portions 
of any trajectory are those in which hydrodynamic interactions with the obstacles are small 
and particles move close to a straight line. As a result, small values of the excess length, 
which correspond to trajectories closer to a straight line, also indicate faster motion of the 
particles, and vice versa. 

D. Suspensions 

Given that the trajectories of a single isolated particle in an obstacle array are not collinear 
with an applied force and may exhibit locking behavior, a natural question is whether this 
behavior persists in the presence of other mobile particles. The basic argument for locking is 
based on the behavior of a single particle approaching a single obstacle but multiple particles 
in general introduce additional interactions which may alter the result. We have therefore 
repeated the above calculations for suspensions of particles. Two cases are considered, both 
using obstacles of radius 2 spaced by 17.2, a "small" simulation involving a 3x3 unit cell 
containing 20 mobile particles of radius 2 and 20 of radius 3, and a "large" simulation with a 
5x5 unit cell, 30 particle of radius 2, 30 of radius 3 and 10 of radius 4. The suspensions are 
fairly dilute, with particle volume fractions of 6.6% and 3.85% for the small and big cases, 
respectively. In Fig. [12] we show the trajectory angle for each size of particle in the two 
systems, when a force of magnitude F = 100 is applied. Although the obstacles evidently 
have an influence on the trajectories, the effect is considerably reduced in comparison to 
the simulations involving a single mobile particle. The interactions between the particles 
alter the motion around the obstacles sufficiently to prevent the accumulation of identical 
individual displacements. On the positive side, there is a clear sensitivity to the particle 
radius in the "big" case, at lower total concentration, which suggest that using flow through 
obstacles as a separation technique is viable for sufficiently dilute suspensions. 
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FIG. 12: (Color online) Trajectory angle vs. forcing angle for suspensions. Left: "small" and right: 
"big" suspension, respectively; with radius 2 (red,+0), radius 3 (green,*) and radius 4 (cyan,o) 
particles displayed separately. 

IV. CONCLUSIONS 

We have used molecular dynamics simulations to expand the range of previous calcula- 
tions of locking behavior in flows past obstacle arrays. The new ingredients of the problem 
considered here are realistic surface roughness, variable particle and obstacle sizes, and a 
comparison of applying a force directly to the particle with forcing the solvent fluid. 

We first considered the behavior of trajectories as a function of Peclet number, and 
observed that even at Pe = 10 the paths of driven particles resembled those at high Pe. At 
general forcing angles, trajectories consisted of segments of variable length aligned along the 
commensurate locking directions of the obstacle lattice, but near the more robust locking 
directions (45° in particular) straighter sinusoidal paths predominated. We also verified the 
theoretical expectation that cylindrical obstacles are needed to produce locking behavior, 
by simulating flow through a three-dimensional lattice of spherical obstacles. We then 
considered the global particle motion for various forcing strengths and orientations and 
for two types for forcing - driving the particle directly and driving the suspending fluid. 
The degree of locking was found to increase with Pe and to be different for the two forcing 
methods, but for Pe = 0(100) did not show the devil's staircase plateaus seen in the Pe — * oo 
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limit. These results were compared in detail to those of Stokesian dynamics calculations. 
Examination of the behavior of particles of different size in the same obstacle lattice suggests 
that smaller particles may lock more effectively. A follow-up study of the motion of multi- 
particle suspensions through obstacles verified that the trajectory angles depend on particle 
size and implies that these may be used as a separation device. Lastly, we considered larger 
particles and obstacles made of more atoms, which leaves the absolute scale of the roughness 
(roughly one atom) unchanged but reduces its scale relative to the particle size, and found 
no significant variation. 

The general result is that independent of these variations, the simulations all show tra- 
jectory locking, in the sense that the average orientation of the trajectory differs from the 
direction in which the force is applied. Clear effects are seen for all sizes and both forc- 
ing mechanisms, once the Peclet number is O(10) or more. However, the sharp step-wise 
transitions which occur in the absence of Brownian motion are smoothed by molecular dif- 
fusion and most of the locking angles observed in the convective limit, shown for example 
in Fig. El are not apparent at finite Peclet number. These results are subject to the usual 
limitation of molecular dynamics calculations: small sizes and short times in comparison 
to many laboratory studies and to continuum methods which do not track atomic motions, 
but the simulations nonetheless exhibit Navier-Stokes behavior and provide relevant results. 
Furthermore, for actual nanoparticle applications the simulations operate directly at the 
relevant scales. 

This material is partially based upon work supported by the National Science Foundation 
under Grant No. CBET-0731032. 
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